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We study two lattice models, the honeycomb lattice (HCL) and a special square lattice (SQL), 
both reducing to the Dirac equation in the continuum limit. In the presence of disorder (gaussian 
potential disorder and random vector potential), we investigate the behaviour of the density of 
states (DOS) numerically and analytically. While an upper bound can be derived for the DOS on 
the SQL at the Dirac point, which is also confirmed by numerical calculations, no such upper limit 
exists on the HCL in the presence of random vector potential. A careful investigation of the lowest 
eigenvalues indeed indicate, that the DOS can possibly be divergent at the Dirac point on the HCL. 
In spite of sharing a common continuum limit, these lattice models exhibit different behaviour. 

PACS numbers: 73.23-b,73.63-b,72.10.Fk 



I. INTRODUCTION 



Graphene, a two-dimensional sheet of carbon atoms forming a honeycomb lattice, has set the stage for studying 
Dirac-type quasiparticles in two dimensional materials^! 3 -. A substantial part of the investigation has been devoted 
to the unusual transport properties of graphene. More recently, also local properties have been studied^. 

Many physical properties depend directly or indirectly on the density of (quantum) states at the Fermi energy. 
Therefore, the density of states (DOS), especially near the Fermi level, is an interesting and important quantity to 
study. Local probing of graphene, such as in the recent STM experiments^, have also raised interest in the local 
DOS. Moreover, the DOS at the Dirac point also plays an important role as an indicator for spontaneous symmetry 
breaking, which causes long-range correlations in graphene^. 

In pure graphene (or for pure Dirac fcrmions), in contrast to disordered graphene, the DOS vanishes linearly like 
p(E) ~ \E\ at the Dirac point E = 0. Scattering by disorder may create new states at any energy, also at E = 0. As 
a consequence, the linear behavior of the DOS at low energies is affected by disorder. On the other hand, the linear 
behavior of the DOS can be considered as a power law of a critical phenomenon with exponent 1. In fact, the phase 
transition in the 2D Ising model is directly linked to this linear behavior of the DOS of 2D Dirac fermions 7 -. A common 
belief is that disorder or additional interaction effects do not destroy the critical phenomenon but only modify the 
exponent of the corresponding power law. This possibility has also been discussed for the Dirac fermions, for instance, 
in the case of a random vector potential^ ' 10 ' 11 . Another possibility is that disorder creates a new intermediate phase 
between the two phases of the pure system^. 

For weak disorder we can apply a perturbation theory with respect to a random vector potential. This approach 
gives a power law 

{p{E))~\E\ a (a<l), (1) 

where the exponent decreases with increasing variance of the disorder distribution g as 

a ~ 1 - g/ir . (2) 

On the other hand, there has been a long debate in the literature whether or not the exponent can have negative 
values for strong disorder (i.e., whether or not there is a divergent average DOS in the case of strong disorder) for the 
model with a single Dirac cone^ ' 11 ' 12 ' 13 . 

The case of two Dirac cones with intervalley scattering has also been discussed intensively in the literatu re 11 ' 15 : 21 ' 22 . 
Intervalley scattering may affect the density of states strongly, leading to a power law with a universal exponent 
a = 1/7 for any strength of disorder—. 

The power law of the density of states has direct implications for the transport properties. The Einstein relation 
states that the conductivity a and the DOS arc proportional to each other: 



a oc p(E)D(E) , 
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where D(E) is the diffusion coefficient. If p(E) vanishes at the Dirac point E = for a > 0, the conductivity also 
vanishes, as long as D(E = 0) is finite. The latter should be the case in the presence of disorder because D(E) 
measures the amount of scattering, since D is proportional to the scattering time r. An exceptional case is a pure 
system, where transport is ballistic (D(E — > 0) — > oo). On the other hand, if p(E) diverges at the Dirac point for 
a < 0, the conductivity also diverges, unless the diffusion coefficient vanishes. 

An alternative approach for the density of states is the self-consistent non-crossing (or Born) approximatio n 16 ! 17 ' 18 . 
The pcrturbativc result of the DOS in Eqs. (fTJ), |(2J) was confirmed for the tight-binding model on the honeycomb 
lattice within the self-consistent calculation^. However, very close to E = an interception of the power law was 
found, indicating a non-zero DOS at E = 0. Moreover, the calculation gave only positive exponents a, even for strong 
disorder, in contrast to the exponent suggested in Ref&ii 

1 - g/n 
a = — . 

1 + g/n 

In order to shed some light on the behavior of the average DOS near the Dirac point, we shall focus in this paper 
on two cases: (i) a single Dirac cone with random vector potential and (ii) the honeycomb lattice with unidirectional 
random bonds. By comparing these two cases we estimate the effect of intervalley scattering on the DOS. 

The paper is organized as follows: After a brief introduction of the tight-binding model for graphene and the 
projection to a single Dirac cone we discuss the underlying symmetries of the models in Sect. 2. Based on these 
considerations we derive a simple expression for the local DOS in the case of the single Dirac cone with random vector 
potential in Sect. 2.2. This allows us in Sect. 2.3 to give an upper bound for the average local DOS. In the second 
part of the paper (Sect. 3) we apply exact diagonalization to the single Dirac cone with random vector potential and 
to the tight-binding model on the honeycomb lattice with unidirectional bond disorder to study the energy levels near 
the Dirac point for finite systems. 

II. MODELS AND SYMMETRIES 

Starting point is a tight-binding model for quasiparticles on the honeycomb lattice. The honeycomb lattice is a 
bipartite lattice. After dividing it into sublattice A and B, the quasiparticles are pseudospin-1/2 particles with respect 
to the two sublattices, and the corresponding Hamiltonian has a chiral symmetry. This allows us to write 

n ~~ 2^ n r,r,r',j' l 'r,j L r',3' > 

r,r' j,j'=l,2 

where r runs over sublattice A and j refers to sublattice A (j = 1) and sublattice B (j = 2). The only energy scale 
of this Hamiltonian is the hopping energy t. Then the Hamiltonian matrix can be expressed with Pauli matrices as^ 

H WL = hiai + h 2 a 2 . (3) 

hi and h 2 , defined on sublattice A, are symmetric and antisymmetric matrices (fv[ = h\, = —h 2 ) 1 respectively. 
The off-diagonal element of the Pauli matrices connect the two sublattices. These properties imply a real symmetric 
Hamiltonian. The corresponding quasiparticle dispersion has two Dirac cones (two "valleys") at low energies. 

A. Dirac Hamiltonian 

Considering quasiparticles at low energies only, we can expand the Hamiltonian around both Dirac points. Then 
we get a model that describes two separate spin-1/2 Dirac spinors. Scattering by disorder can, in principle, connect 
these two Dirac cones (valleys). It has been discussed that this leads to the SU (2) Wess-Zumino-Witten models (but 
see also RefJ^). On the other hand, if inter-cone scattering is ignored (for instance, by assuming a smooth scattering 
potential that is constant on the scale of the lattice spacing), the two valleys of the model are completely isolated 
from each other and each valley can be studied separately. Then disorder can appear as a random scalar potential, a 
random mass or a random vector potential 8 -. Only the latter preserves the continuous chiral symmetry. It is believed 
that this type of disorder is related to ripples in the graphene shee t 19 ' 20 . The corresponding Hamiltonian Ho is again 
a chiral spinor-1/2 Hamiltonian but, in contrast to the real symmetric tight-binding Hamiltonian on the honeycomb 
lattice -ffHCL, it breaks the time-reversal invariancc 



H D = h\(Ji + h 2 a 2 



(4) 
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hi and /12 are now antisymmetric spatial matrices (hj — —hj (j = 1, 2)) with imaginary matrix elements, and a here 
denotes the physical spin, this is why this Hamiltonian breaks the time reversal invariance. This gives H D = a^H-DCi- 
Moreover, we assume that hj are lattice hopping matrix elements with nearest-neighbor elements on a square lattice 
whose continuum limit is the j component of the 2D gradient Vj. This fictitious square lattice is sketched in Fig. [1] 
with spin dependent hopping amplitudes. Thus, the Hamiltonian Hd describes lattice Dirac fermions. The lattice 
constant is not that of the original honeycomb lattice of the graphene sheet but larger, and related to the projection 
onto a single Dirac cone. In this respect the lattice structure of Hd corresponds with the network approximation of 
the honeycomb lattice^. 



FIG. 1: (Color online) The square lattice, whose continuum limit is the Dirac Hamiltonian is visualized. Filled red and empty 
black circles denote up and down spins at a given lattice point, thick/thin lines denote the hopping/lattice. The hopping matrix 
elements are indicated. Note the spin dependent hopping amplitudes! 

In the following, disorder due to ripples will be considered. This can be represented by a random vector potential 
(Vi,r, V 2 , r ) as 



H ={h 1 + Vi)(Ti + O2 + V 2 )a 2 



(5) 



This Hamiltonian has three essential symmetry properties: It is Hcrmitian (i.e. W = H), and it satisfies the following 
relations: 



a 3 Ha 3 



-H 



(6) 



and with the staggered diagonal matrix D 



(-lr+^.^v 



we get (cf. Appendix A) 

aiDH T Dai = H . 

The fact that H is Hcrmitian implies for the Green's function G(ie) = (ie + H)^ 1 the relation 

G\ie) = G(-ie) . 

Moreover, Eq. © implies 

a 3 G(ie)a 3 = — G{— ie) , 

and Eq. ([7]) implies 

a x DG(it) T BGi = G(ie) . 

The spatial diagonal elements of the Green's function G rr (ie) can be expressed in terms of Pauli matrices as 

G rr (ie) = go(ie)a + gi{ie)ai + 92^)02 + 93(^)^3 ■ (11) 
The three relations in Eqs. © - (|10[) provide the following relations between the coefficients of the Pauli matrices: 
9o( ie ) = 9o(~ie) = -ffo(ie) , 3*(* e ) = ffi(-ie) = 3i(« e ) > 92( ie ) = 92(-ie) = g 2 (ie) , 

93^e) = . 

Note, that this is a clear consequence of Eq. (fT0| , which holds true only on the square lattice. Thus, 50 is purely 
imaginary, whereas gi and g 2 are real and (73 vanishes: 



(7) 
(8) 
(9) 
(10) 



G rr (ie) = g (ie)a + gi{ie)ai + g 2 («e)cr 2 . 



(12) 
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B. Local density of states of Dirac fermions 

The Green's function G = (ie + Hjj)^ 1 allows us to write for the local DOS for a fixed random disorder configuration 



Pr 



2tt 



ImTr 2 (G r 



(13) 



where e > is implicitly sent to zero, and the Tr is taken over the Pauli matrices. As a function of the random vector 
potential at site r (Vi, r , Vz, r ), the local DOS p r of the Green's function in Eq. fT2|) has a Lorentzian form (cf. Eq. 
([Bl~]) in Appendix B): 



Pr 



(A + 6) 



tt (A + e) 2 + [X x + V, r ) 2 + (A 2 + V 2 , r f 



(14) 



with some real variables X±, A 2 and a positive real variable e + A , where the latter is proportional to e. They depend 
on V\y , V2. r ' for r' =/= r but not on V\ >r , V 2 _ r . This expression can also be used to determine the DOS away from the 
Dirac point at energy E ^ by replacing e — > e — iE: 



pJE) = -Re 

7T 



(A 



iE) 



(A + e - iE) 2 + (Aj + V hr ) 2 + (A 2 + V 2 , r ¥ 



(15) 



It should be noticed that this form of the local DOS is very special for the Green's function in Eq. (fT2|) . For instance, 
we would not get a Lorentzian in the case of a random scalar potential. 

Expression (fl5|) enables us to evaluate the local DOS p r {E) for an impurity at site r. According to Eq. (|B2[) the 
parameters Xj (J = 0, 1, 2, 3) of the system without disorder are 



A = -e + iE + i 



9 9 : 

.91 - 9o 



Ai — Xo — At — 



where 



30 



ie + E 



d 2 k 



(e - iE) 2 + k 2 (2tt) 2 



(16) 



The local DOS of Eq. JT5J) then reads 



p r (S) = -Re 

TT 



(l + g 1 V 1 , r ) 2 -g 2 V 2 r 



We can also study a local scalar potential E r by adding the latter to the energy E in go of Eq. (fH)|) . The contribution 
of the local potentials E r and Vi >r to (p r (E)) is quite different, as shown in Fig. [2] While the scalar potential creates 
mostly states at and very close to the Dirac point, the vector potential creates states in some distance from the Dirac 
point. 

A direct evaluation of the variables Xj (J = 0, 1, 2, 3) is difficult in the general case, where we have a random vector 
potential at all sites. However, for finite and sufficiently small systems an exact diagonalization is possible. Moreover, 
we can derive an upper bound for the average local DOS. This will be discussed in the next section. 



C. Upper bound for the DOS of Dirac fermions 

Now we perform the integration with respect to (Vi , V 2 ) for all sites to evaluate the average local DOS. For simplicity, 
we consider only the Dirac point E = here: 



(Pr}= J PrHP(V iy )dV h ,,P(V 2 , rl )dV 2 , rl 



(17) 



First, we perform the integration with respect to Vi, r , using the expression of p r in Eq. (|14[) 

(An + e) z + (Ax + Vi >r y + (A 2 + v 2 ,rr 




lo/D 

FIG. 2: (Color online) Average local DOS (p r (E)) of the Dirac Hamiltonian (Eq. @) for a local random vector potential Vi, r 
(red curve) and a local random scalar potential E r (blue dashed curve). The potentials are box distributed with — 1 < Vi. r < 1 
and —0.1 < E r < 0.1, D is the cutoff in the continuum theory. 

An upper bound for this integral is obtained from pulling out the maximum of the distribution density P{V\. r ) which 
we call P m : P(Vi >r ) < P m . This gives 



J p r P(Vl, r )dV lt r < ^ J 



(Xo + e) 



7dV hr , 



+ e) 2 + {X l + Vi. r ) 2 + (X 2 + V 2 , r f 

and after integrating over the Lorentzian function, which gives tt, the right-hand side becomes P„ 

J PrP(Vi,r)dVl,r < Pm ■ 
Going back to the expression in Eq. (fT"7|) . we obtain 

{pr) < P m J P{V 2 ,r)dV 2 ,r J Y[ P(V L , rl )dVl,r> P(V 2 y)dV 2 y = P m . 

In other words, the averaged local DOS at the Dirac point E = has an upper bound: 

(Pr) = ^-Tr 2 ((lmG rr )) < max P(V) . 

Z7T — oo<\/<oo 



(18) 



This means that for any smooth bounded distribution of V\ iT (e.g. for a Gaussian) the corresponding average local 
DOS p r is finite. For discrete distributions, such as a binary alloy, the upper bound is infinite though. 

III. EXACT DIAGONALIZATION 

For a better understanding of the details of the DOS, we employ an exact diagonalization study on small clusters 
for both models, the Hamiltonian of Eq. §5§ on the original honeycomb lattice (HCL) and the Hamiltonian of Eq. ([5]) 
on the effective square lattice (SQL). Although both models reduce to the same continuum limit of Dirac fcrmions 
with random vector potential, they possess distinct structures in the DOS, as we will discuss below. We use Gaussian 
disorder with standard deviation V (i.e., V 2 is the variance). 



1. Density of states by ED 

Determining the DOS of the infinite system by studying a finite system is a difficult task, since any finite system 
possesses distinct energy levels, resulting in separate Dirac delta peaks in the DOS at the quasiparticle energies. The 
DOS becomes continuous only in the thermodynamic limit. In order to avoid this problem, we choose an indirect 
approach to evaluate the DOS by counting the number of eigenvalues in a narrow frequency range around a given 
energy E. Strictly speaking, this leads to the number of states around E, but if the DOS is a smooth function, this 
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provides us with a sensible definition. We obtain the DOS shown in Fig. [3] on a 100 x 100 HCL cluster with periodic 
boundary conditions for unidirectional bond and potential disorder, using a t/500 wide energy windows, where t is the 
uniform hopping amplitude. For comparison, we also show the result of the self-consistent non-crossing approximation 
(SCNCA) on the HCL^S. As is seen, the agreement is surprisingly good for weak disorder, except for the case of 
bond disorder in a very close vicinity of the Dirac point. There, for V\ < 0.6i, the residual DOS remains zero, which 
is in contrast to the finite, although exponentially small, residual value for the case of potential disorder, described 
correctly by the SCNCA. A narrow peak appears at the Dirac point (DP) for bond disorder if V\ > 0.6i. Whether this 
peak remains finite or diverges cannot be decided within this calculation of the DOS. It should be mentioned that the 
DOS on a SQL is qualitatively similar to the potential disorder case on a HCL for strong disorder. In particular, it 
never diverges at the DP. The anomalous behavior close to the DP is obvious in perturbation theory as welli^, where 
a dynamically generated low energy scale, similar to the Kondo scale, separates the high and low energy regions in 
the DOS. 




1 2 1 2 1 2 

Lojt Lu/t U)/t 



FIG. 3: (Color online) The DOS is shown as obtained by exact diagonalization on 100 x 100 honeycomb clusters with Gaussian 
unidirectional bond disorder (left panel), potential disorder (middle panel) after 1000 averages for Vi/t = 0.3 (blue), 0.5 (red), 
0.7 (black), 0.9 (magenta) and 1 (green). The right panel shows the corresponding self-consistent non-crossing approximation 
for the same parameters for the HCL. The inset shows the narrow peak at the DP for the unidirectional case. The SCNCA 
leads to the same result for pure unidirectional bond or potential disorder. Note the nice agreement between the numerical and 
analytical results for weak disorder! 



2. Eigenvalues 

The investigation of the lowest eigenvalues in the case of unidirectional bond disorder, determining the residual 
DOS, may reveal some structures which arc responsible for the aforementioned behavior of the DOS near the DP. 
Therefore, we take a single disorder realization of Hdi S , chosen randomly according to a Gaussian distribution. Then we 
diagonalizc HiiCL + ViHdis, using the Lanczos algorithm, and retain the 200 eigenvalues closest to the DP (symmetric 
to the DP). This procedure is repeated for different values of V\. The result is shown in Figs. [4] and [5] as a function 
of the disorder strength for a 1000 x 1000 cluster on the HCL and a 708 x 708 cluster on the SQL, having almost 
exactly the same number of states. This reveals three different regimes: 

(i) for weak disorder, the distribution of the eigenvalues is rather dilute and is not influenced significantly by 
disorder. This can explain the zero residual DOS in this case, where a slight rearrangement of the eigenvalues change 
only the slope of the vanishing DOS. 

(ii) Around V\ ~ 0.7i, the pattern changes drastically for the HCL, where the spectrum becomes very dense close 
to zero energy. It keeps on decreasing monotonically down to zero energy. This behavior is responsible for the peak 
and a possible divergence of the DOS. 

(iii) For strong disorder {Vi/t ~ 5), the eigenvalues depart from the DP again. This crossover is related to finite 
size effects, since the characteristic disorder value shifts markedly to higher values with increasing system size. This 
is different for the SQL. At low values of Vy, the DOS behaves similarly for the HCL as well as for the SQL, where the 
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DOS goes down in a power-law fashion, with decreasing exponent, but retains a finite value at the DP. For V\ > t, 
however, the eigenvalue pattern is strongly affected only on the HCL by the explicit value of the disorder. A direct 
study of the DOS reveals no peak around the DP for the SQL but a finite residual value. This reflects the upper 
bound which was derived in Sect. 2.3. 

In order to obtain the DOS, we employ another approach for evaluating this quantity at the DP, which was 
introduced in Ref. Hfl We determine the number of states N(E) in a given energy interval E around the Dirac point 
and define the DOS as lim^^o N(E)/E. As is seen in Fig. \7\ the resulting DOS for Vi > Q.7t shows an upturn 
with decreasing energy for bond disorder, which may be indicative for a diverging nature of the DOS. The DOS for 
V\ = 0.5t still goes to zero, but the 0.7 data increases monotonically with decreasing energy. This supports the picture, 
that the residual DOS is indeed zero for V < 0.6...0.7i, and changes to a diverging behavior afterwards. The results 
for V = 0.3t are probably strongly affected by finite size effects. By fitting the resulting curves with a power law, we 
determine the exponents (a) which is characterizing the DOS close to the DP (cf. Fig. [5]). From a the dynamical 
exponent z follows as z = 2/(1 + a). According to Refill, the latter changes its behavior at z = 3, which is reached 
here at Vi/t ~ 2.5, and it increases linearly with Vi. For comparison, the case of potential disorder is plotted as well, 
where the DOS tends smoothly to a constant value at E = 0. The SQL with V\ disorder exhibits qualitatively similar 
behavior to the potential disorder case on the HCL. 




1 2 3 4 5 6 

Vi/t 



FIG. 4: (Color online) The evolution of the lowest 100 eigenvalues above the DP is shown for a 1000 x 1000 HCL cluster with a 
given Gaussian disorder configuration on a semilogarithmic scale, by changing the strength of the disorder. The inset enlarges 
the low energy structures and the transition from vanishing to diverging behavior. For Vi > 0.7T, the eigenvalues start to 
approach zero rapidly, as is obvious from the semilogarithmic scale. Their increasing behaviour for Vi > 5 is due to finite size 
effects. The statistics of the eigenvalues at Vi = 3t is depicted in Fig. [6] 



3. Finite DOS on the SQL 

Now we turn our attention to the square lattice model in Eq. (5). For the pure system, there is no difference 
between the HCL and the SQL for the DOS near the DP, since excitations close to half filling are Dirac fermions 
in both cases. Thus, the DOS increases linearly with energy. It also exhibits a weak logarithmic singularity at the 
saddle point of the spectrum, and falls off monotonically with increasing energy towards the band edge, as is seen in 
the inset of Fig. |TD] The V\ disorder in Eq. [5] on the lattice model plays the role of a random vector potential, which 
is perpendicular to the (pseudo) spin quantization axis 173. In the presence of V\ disorder the DOS on the SQL is 
different from that of the HCL with unidirectional bond disorder: no peak develops at zero energy for strong disorder, 
and the DOS terminates at a finite value with vanishing slope, similarly to potential disorder in the HCL. Using an 
energy window of t/500 as for the HCL, we can evaluate the DOS as described above. The residual values are plotted 
in Fig. 1101 and compared with the upper bound. As is seen, the upper bound becomes very sharp for strong disorder 
in this case, and does not seem to apply to the HCL with a possibly diverging DOS. 
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FIG. 5: (Color online) The evolution of the lowest 100 eigenvalues above the DP is shown for a 708 x 708 SQL cluster with a 
given Gaussian disorder configuration on a semilogarithmic scale, by changing the strength of the disorder. The inset enlarges 
the low energy structures and the transition from vanishing to diverging behavior. As opposed the the HCL, the structure of 
the eigenvalues hardly changes for Vi > t. Their distribution is shown in Fig. [6] 




1 2 3 4 5 6 



Ei *10- 5 

FIG. 6: (Color online) A typical distribution of the lowest eigenvalues is shown for V/t = 3 for both, the HCL and the SQL. 
In the former case, the eigenvalues precipitate to zero very fast, resulting in a sharp peak around zero energy. As opposed to 
this, the distribution for the SQL is more uniform, yielding a nondiverging constant DOS. 

IV. CONCLUSIONS 

We have evaluated the eigenvalues and the average DOS for the tight-binding model on the honeycomb lattice with 
random unidirectional bonds and for Dirac fermions on the square lattice with random vector potential. Both models 
have the same continuum limit, namely Dirac fermions with a random vector potential. However, in their lattice form 
they differ substantially near the Dirac point: In the model on the honeycomb lattice the average DOS has a sharp 
peak which is not present in the model on the square lattice. Although it is not entirely clear, whether or not this 
peak survives the limit of the infinite system, its existence on the finite cluster is remarkable. The evolution of the 
eigenvalues close to the Dirac point in large systems supports the idea of a diverging peak in the DOS. 

We have studied the effect of potential disorder on the honeycomb lattice as well, which exhibits qualitatively similar 
behaviour to the square lattice with random vector potential, but differs from the case of random unidirectional bond 
on the honeycomb lattice at low energies in the DOS: the residual DOS always takes a finite, although exponentially 
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FIG. 7: (Color online) The number of states divided by energy (~ p(^)) is plotted as obtained by exact diagonalization on 
100x100 honeycomb clusters with Gaussian unidirectional bond disorder (left panel) and potential disorder (right panel) after 
1000 averages for several values of the disorder. The upturn with decreasing energy for bond disorder is indicative to the 
diverging DOS at E = for V/t > 0.6. 




FIG. 8: (Color online) The number of states divided by energy (~ p(^)) is plotted as obtained by exact diagonalization on 
90x90 square lattice after 1000 averages for several values of the Vi disorder. It resembles closely to the potential disorder case 
of the HCL. 

small value. These results can surprisingly well be reproduced for weak and moderately strong disorder using the 
self-consistent non-crossing approximation, expect for the low energy structures in the case of bond disorder. 

Using the mapping of the model on the HCL to the SU(2) gauge field theoi y 11 ' 15 ^ 1 ' 22 , the presumably exact power 
law of the latter p ~ of Refpii represents a puzzle for the approximation of disordered lattice models by their 

corresponding continuum counterparts. The same is true for the model on the square lattice, where the DOS at the 
Dirac point has an upper bound according to Eq. (|18[) . In contrast, for the continuum limit several groups found a 
power law with the exponent&ii 

1 - g/rr 

a = ~i T 

1 + g/ir 

which is negative for sufficiently strong disorder. This poses severe questions on the applicability of universality idea. 
Although both lattice models reduce to the same continuum limit and arc expected to behave in a similar manner, 
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Vi/t 

FIG. 9: (Color online) The exponents of the DOS (p(w) ~ <^ Q ) and the dynamical exponent z = 2/(1 + a) axe plotted for the 
HCL for strong disorder. Note the horizontal axis, which is the standard deviation and not the variance. 
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Vi/t 



FIG. 10: (Color online) The residual DOS of the square lattice, Eq. [5] is plotted obtained on 90 x 90 and 30 x 30 clusters with 
Gaussian Vi disorder (red squares and blue circles) after 10 3 and 10 4 averages, respectively. The black straight line is the upper 
bound, p(0) < l/\/27rVi. Note, that it does not involve any fitting parameter, and becomes very sharp for strong disorder. 
Larger systems show similar behavior. Inset: the DOS of a 90x90 SQL is shown after 1000 averages for Vi/t—Q.5, 0.8, 1, 1.5, 2 
and 2.5 with decreasing peak-position at uj = 2t. 



as dictated by the common continuum limit, this is apparently not the case here. We have also checked the case of 
uniform disorder distribution, and found similar results. The above results were found to be robust with respect to 
variations of system size, boundary conditions, and disorder distribution. 
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APPENDIX A: DISCRETE SYMMETRY 

From hj = ~hj and a\ — a\, <r 2 = — a 2 follows 

H T = (-h X + V^CT! - (-h 2 + V 2 )ot 2 . 
Next, D changes the sign of nearest-neighbor matrix elements: 

DH T D = (ht + Vifa - (h 2 + V 2 )cj 2 , 

and o\ anticommutes with a 2 : 

a x BH T Da x = (h x + V l )(J l + (h 2 + V 2 )a 2 = H . 



APPENDIX B: MATRIX ELEMENTS OF THE GREEN'S FUNCTION 

The spatial diagonal matrix elements of the Green's function have been given in Eq. (|12p . Another way to write 
G rr is by projecting it with P r onto the site r. This gives the matrix identity^ 



P r GP r = [it + V r at + V;<r 2 - P r H{l - P r )Gi_ Pr (l - P r )HP r ]p l 



where Gi_p r is the Green's function G(ie) = (ie + H) 1 on the Hilbert space where the site r has been removed. The 
2x2 matrix P r H(l — P r )Gi-p r (l — P r )HP r does not depend on the random variables V r and VJ. Its general form is 



P r H(l - P^d-pjl - P r )HP r 



Therefore, G rr reads 



iX Q + X 3 —iX 2 + Xi 
iX 2 + X\ iX a — Xs 



Gv 



it + iX Q + X 3 
lX 2 + Xt + V r + iV' 



-lX 2 + X 1 + V r - %V{ 

it + iX - X 3 



1 



(e + X ) 2 + XI + {X l + V r ) 2 + (X 2 + V r r f 
This result can be compared with Eq. (TL^)) to obtain the relations 

9i v 9o 



it + iX - X 3 iX 2 + Xi + V r + iVi 
iX 2 + X l + V r - iV r ' it + iX + X 3 



X X = -V r + 



9o + 9l + 92 



iXn = —it — 



-So + 9i + 92 



X 2 = -VI 



.92 



-9o +9i+92 ' 



(Bl) 



(B2) 



12 



and 

X 3 = . 

All three matrix elements Xq, X\, X 2 are real, since go is purely imaginary, and gx as well as g2 are real. 

Finally, we can use the block-matrix inverse to show that go is proportional to — it with a positive proportionality 
factor. Choosing the diagonal blocks with respect to the sublatticc (or spinor) index j, we obtain 

Gn = [it - (hi +V X - ih 2 - iV 2 ){h l +V X + ih 2 + iV 2 )/it)]- 1 = -it[t 2 + {hx + Vx- ih 2 - iV 2 )(hx + V 1 + ih 2 + iV 2 )]^ 

= -te[e 2 + (hx +Vi- ih 2 - iV 2 ){hi +Vi- ih 2 - ^V 2 ) t ^ 1 

and 

G 22 = [it - (h x + Vx + ih 2 + iV 2 )(hi +V X - ih 2 - iV^/it)]- 1 = -it[t 2 + (hx + Vi + ih 2 + iV 2 )(h x +Vi- ih 2 - iV^]- 1 

= -it[t 2 + (hx +Vi- %h 2 - W 2 f(h x + Vx - ih 2 - iVz)}- 1 . 
Thus iX + it = cit with c > 0. 



